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Abstract. We provide a convergence analysis for a new fractional time-stepping technique for 
the incompressible Navier-Stokes equations based on direction splitting. This new technique is of 
linear complexity, unconditionally stable and convergent, and suitable for massive parallelization. 



1. Introduction 

This work is concerned with the analysis of a new class of approximation techniques for the 
solution of the time-dependent incompressible Navier-Stokes equations based on direction splitting. 
This new technique requires, independently of the space dimension, only the solution of a sequence 
of one-dimensional problems, thus having linear complexity. The main claims of this paper are 
that this technique is unconditionally stable and supcrlincarly convergent with respect to the time 
discretization parameter and is suitable for massive parallelization. 

We consider the Stokes equations written in terms of velocity u and pressure p on a finite time 
interval [0, T] and in a cubic domain fi = (0, l) d with d — 2 or 3: 



(1.1) 



'u t -Au + Vp = /, infix(0,T], 

Vu = 0, infix [0,7*1, 

u\an = 0, in (0,T], 

u| t= o = u , in fi. 



where / is a smooth source term and Uo is a solenoidal initial velocity field with zero normal trace. 
The nonlinear term in the momentum equation of the Navier-Stokes equations is not accounted for 
since it does not interfere with the incompressibility constraint. The fluid density is assumed to be 
constant and has been put into the normalization constants. 



Once time is discretized, (1.1) reduces to a generalized Stokes system at each time step. Solving 
this coupled system often proves computer intensive and is not easy to solve efficiently in parallel 
due to the saddle point structure induced by the incompressibility constraint. Alternative more 
efficient approaches consist of uncoupling the velocity and the pressure using so-called projection 
algorithms. 
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Projection algorithms date back to the late 1960s and stem from the seminal works of Chorin 
[2] and Temam [23]. These methods and various improvements thereof are still, to the best of 
our knowledge, the methods of choice in the CFD community. Although in the 1980s and 1990s 
these techniques underwent some evolution and their properties are now fairly well understood 
[T51 [201 122 1221 12S1 EH H] (the reader is referred to [5] for an overview), the same fundamental idea 
of decomposing vector fields into a divergence-free part and a gradient has remained unchanged 
over the years and has been challenged only recently in [12]. For all these schemes, the total cost 
per time step is that of solving one vector- valued advection-diffusion equation and one scalar- valued 
Poisson equation with homogeneous Neumann boundary conditions. For very large size problems, 
the cost of solving the Poisson equation is dominant. To address this issue, Guermond and Minev 
have proposed a new method in [5]. The main idea consists of abandoning the projection paradigm, 
as in [T23 j an d replacing the Poisson equation by a direction splitting strategy. This requires to solve 
a sequence of one-dimensional elliptic problems instead of one multidimensional Poisson equation. 
The first-order accurate variant of method has been shown to be unconditionally stable in [9] . 

In this paper we pursue further the ideas introduced/announced in [9] in the sense that in addition 
to splitting the pressure-correction, we also apply a direction splitting technique to the momentum 
equation, thus further reducing the overall computational cost of the method. We prove that the 
totally split method is convergent and we provide error estimates. 

Applying direction splitting to the momentum equation is not a new idea. For instance, in |24[ 
Section 3.7.2] Temam studies a projection method where the solution of the momentum equation 
is obtained using direction splitting and the incompressibility constraint is enforced by means of a 
Poisson equation. Stability and convergence of the scheme are proved therein but no error estimates 
are provided. Lu, Neittaanmaki and Tai show in JTHl \H\ that this scheme is C( r5 ) accurate, r 
being the time-step. Our work differs from these previous results mainly in two directions. First, 
we adopt a direction splitting strategy for the computation of the pressure-correction which renders 
the method extremely fast and massively parallelizable. Second, we provide error estimates for the 
proposed scheme, and we show that the so-called standard version of the scheme is C(r)-accurate 
in all quantities irrespective of the space dimension and the rotational version is 0{t s )-accurate in 
two space dimensions. Numerical experiments show that the result holds true also in three space 
dimensions and the actual convergence rate is higher than £)(r 5 ) in two and three space dimensions. 
The algorithm has been implemented in a parallel code which has been observed to have optimal 
weak scalability. This code has been used to compute the transient regime on the three-dimensional 
lid-driven cavity at R e = 1000 and R e = 5000 on a mesh composed of 2 10 9 nodes on 512 processors 
only. 

This paper is organized as follows. Section |1.1| introduces the notation and establishes some 
preliminary results. The new algorithm is described in Section [2| two-dimensional and three- 
dimensional variants of the algorithm are presented in §2.1| and §2.2[ respectively. The convergence 
analysis of the standard form of the algorithm is done in Section [3] and the analysis of the rotational 
form is done in Section [4] In Section [5] we briefly discuss the BDF2 technique to march in time. 
Finally, we present numerical experiments in Section [6] to illustrate the performance of this new 
class of algorithms. 



1.1. Notation and Preliminaries. We consider the time-dependent Stokes system (1.1) on the 
finite time interval [0,T] and in the cubic domain ft := (0, l) d with d — 2 or 3. 

We henceforth consider only the time discretization of the system to simplify the discussion. 
Handling the space discretization is a secondary issue, and the reader is referred to [7[ QT] for the 
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techniques that can be used for this purpose. Let r > be a time step (for simplicity taken uniform) 
and let tk = kr for < k < K := \T/t\. Let E be a normed space, with norm || • \\e- For any 
time-dependent function ip : [0,T] — > E, we denote tp k :— ip(tk) and the sequence {ip k }k=o,...,K is 
denoted by ip T . To simplify the notation we define the time-increment operator 5 by setting 

(1.2) 5ip k : =V fe -V> fc ~\ 
and the time- average by 

(1.3) ^ := 
We also define the following discrete norms: 
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(1.4) kEll^HM ' Il^r||/« (B) := n^c {||^ fc || B }. 

V fe=o / 

The space of functions tp : [0,T] — > E that are such that the map (0,T) 3 t — > HV'WIIe e M is 
LP-integrable is indifferently denoted L p ((0, T);E) or U>{E). 

No notational distinction is done between scalar or vector-valued functions but spaces of vector- 
valued functions are identified with bold fonts. We use the standard Sobolev spaces W m ' p (d), for 
< m < oo and 1 < p < oo. The closure with respect to the norm || • ||vK m >p of the space of 
C°°-functions compactly supported in f2 is denoted Wq™' p (Q). To simplify the notation, the Hilbert 
space W s ' 2 (fl) (resp. W^ 2 (rt)) is denoted H'(Q) (rcsp. i? s (O)). We define Lf =0 (Q) (rcsp. H]^(Q.)) 
the space that is composed of those functions in L 2 (Q) (resp. H 1 ^})) that are of zero mean. The 
scalar product of L 2 (f2) and Lj =0 (Q) is denoted (•, •) and we define 

(1-5) |M|*i := ||Vg|| L2 , Vq€Hl =0 (Q), 

(1.6) ||v|| H i := ||V^|| L 2, WeHj(Q). 
Finally we recall that 

(1.7) = ||Vt;||| 2 + ||Vx«!| 2 j2 , V.eHjfH). 

Henceforth c denotes a generic constant whose value may change at each occurrence. This 
constant may depend on the data of the problem and its exact solution, but it does not depend on 
the discretization parameter r or the solution of the numerical scheme. 

1.2. Direction Splitting Pressure Operator. We assume that we have at hand an operator 
A : D(A) C Lj =0 (£l) — > Lf =0 (£V) which is unbounded, closed and satisfies 

(1.8) ||Vg|| 2 2 < (Aq,q), VqeD(A). 
This property implies that the map D(A) 3 q \\q\\A € M where 

(1.9) \\q\\A:=(Aq,q)K VqeD(A), 
is a norm. We also define the scalar product 

(1.10) (p,q) A :=(Ap,q), Vp,qeD(A). 

A natural example for A consists of using A = — Ajv, where A^r is the Laplace operator sup- 
plemented with homogeneous Neumann boundary conditions. This operator is the workhorse of 
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classical projection methods. The main originality of the method that we are going to consider 
consists of introducing a direction factorization of this operator. In two space dimension we define 



(1.11) 



(1.12) 



A:= (l-d XX )(l-dyy), 

D(A) := {p e LjLo(fi) : d w p,Ap £ L 2 (Q) : d yP \ y=0 .i = 0, 0„(1 - d yy )p\ x=0 , 
and in three dimensions 

' A:= (l-d xx )(l-d yy )(l-d zz ), 

D{A) := {p e Lf =0 (n) : d zz p, (1 - d yy )(l - d zz )p, Ap e L 2 (ft) : 

d z p\ z= oA = 0, d„(l - d zz )p\ y = ,i = 0, 9 X (1 - d yy ){l - d zz )p\ x=0 ,i = 0} , 

The graph norm is denoted || • ||.d(a) both in two and three space dimensions. 

Proposition 1.1. The operator A defined in ( |1.11[ ) or (1.12), in two or three space dimensions, 
respectively, satisfies (1.8). 

Proof. See [9]. □ 



d x pi\x=o,i 


= 0, 


dyP2\y=0,l 


= 0, 


d z p\ z =o,i 


= 0. 



One interesting feature of the operators defined by ( 1.11 ) and ( 1.12 ) is that solving the equation 
Ap = f for / € L 2 (Q) only requires to solve one-dimensional problems. For instance, the solution 
of Ap = f in three space dimensions is obtained by solving for p±, p2, and p so that 

Pi - d xx pi = /, 

P2 - d yy p 2 = Pi , 

p- d zz p = p 2 , 

Finally we introduce the Hilbert space Y to be the completion of the space of smooth scalar- 
valued functions with respect to the norm || • H^: 

(i.i3) Y-.= e°°(ct) HA nLf =0 (n). 

The extension of the scalar product (•, ■) a to Y is abusively denoted (•, -)a- For instance if A is 
defined as in ( |l.ll[ ) or (1.12), the space Y is characterized as follows: 

flu) y= UqeH^ (n);d xy qeL 2 (n)}, in M 2 , 

\{q£Hl =0 (n);d xy q,d yz q,d zx q,d xyz qeL 2 (n)} inM 3 . 

Note that the boundary conditions associated with D (A) have disappeared from Y and the || • ||^- 
norm (which is also the norm in Y) is characterized by 

\d xy q\\h, 

vyqWh + \\d yz q\\ 2 L 2 + \\d zx q\\ 2 L2 + \\d xyz q\\ 2 L2 , 



(1.15) 




in 

in 



1.3. Direction Splitting Velocity Operator. To be able to handle the two-dimensional and 
three-dimensional error analysis in a unified framework we introduce the following unbounded 
closed operator 



(1.16) 

with domain 
(1.17) 



Bv := 



^xxyyV 
(,$xxyy 



$yyzz ~t~ d z 



2 ^xxyyzz)^ 



D{B) := {v e Hj(n); Bv G L 2 (ft)}. 
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The graph norm is denoted 



\D(B)- 



Lemma 1.1. The bilinear form D(B)xD(B) 3 (v,w) 
the following holds for all v E D(B): 



(v,Bw) E K is symmetric positive and 



(1.18) 



(v, Bv) 



|ftcj/ u |l:L,2 

I^Ht 2 + ||ft, Z u|| 2 2 



\d xz v» 2 



L 2 



511ft 



in R 2 , 
in K 3 . 



Proof. Let us consider the two-dimensional case first. Using the Fubini-Tonelli Theorem and inte- 
grating by parts repeatedly we obtain 



v=o 



v=i 



d y vd yxx v dy 



v=o 



d.r 



y=o 
= \\d xy v\\i2. 



dyVdy X V 



X=l 

x=0 



x=l 



x=0 



{d yx v) 2 dx 



dy 



Note that we used v\ y= o t i = and d y v\ x =Q > i = which is a consequence of w| x =o,i = 0. The 
three-dimensional result is obtained similarly; the details are left to the reader. □ 

To simplify notation we now define the norm 

(1.19) \\v\\ B :=(v,Bv)2, v E D(B), 
and we define the following Hilbert space 

(1.20) Z :=C°°(fi) IMlB nHj(fi). 

The extension of the scalar product (•, -)s to Z is abusively denoted (•, •) b- 

1.4. The Right-Inverse of the Stokes Operator. To describe solenoidal vector fields we intro- 
duce the classical spaces 



(1.21) 



H:={oe L 2 (fl) : V-v = 0, vn\ m = 0} , V:=Hfl H^tt), 



where n is the outer unit normal to d£l and we denote by Ph the L 2 -projection onto H. It is also 
useful to introduce the right-inverse of the Stokes operator S : L 2 (S1) — > V defined as follows: for 
any / e L 2 (f2) we denote (Sf, q) E V x Lf =0 (Q) the pair such that 

-ASf + Vq = f, in O, 

(1.22) { VSf = 0, in O, 

sf = o, on an. 

Given the particular domain that we consider in this work, the inverse Stokes operator is bounded 
from L 2 (£l) to H 2 n V (cf. [3]), i.e., ||5/||h 2 < c||/||l 2 - Finally, we introduce the seminorm 

(1.23) \v\l:=(Sv,v), VvEL 2 (n), 

and, we recall (see e.g. [3H [7J H3]), that for every 7 E (0, 1), there exists c(j) > so that the 
following holds for every v E Hq(£1) 

(1.24) (WSv,Wv) > (l-j)\\v\\h-c( 7 )\\v-P H v\\l 2 . 
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2. Description of the Scheme 

We describe the direction splitting algorithm in two and three space dimensions in this section. 
The stability and convergence analysis is done in the subsequent sections. 

2.1. Two Space Dimensions. To simplify the presentation, we assume for the time being that 
the space dimension is two (d — 2) and we defer to [ ]2.2| the discussion of the three dimensional 
case. 

The scheme computes three sequences of variables {u k }, {cf> k ~^}, and {p k ~?} that approximate 
the velocity, the pressure-correction, and the pressure, respectively. 

• Pressure predictor: Denoting by p the pressure field at t — and </>* ,_ 5 an approximation 
of |tc^p(0), the algorithm is initialized by setting p~^ = p and (f>~i = <fr*'~i . Then, for 
all k > a pressure predictor is computed as follows: 

(2.1) p*' k+1 i =p*-a +0 fc "5. 

• Velocity update: The velocity field is initialized by setting u° = Uo, and for all k > the 
velocity update is computed by solving the following series of one-dimensional problems: 
Find w fc+ 2 and u k+1 such that 

(2.2) u k +^-u k _ d ^ kH _ Q ^ k + Vp±ikH = fkH uk+ , ^ = 0j 

fc+1 fe+" 

(2.3) - d xx u k+ i d vy u k+1 + Vp*' k+l > = f k+ i , « fe+ Vo,i - 0. 

• Penalty step: The pressure-correction fe+ 2 is computed by solving 

(2.4) A4> k+ ^ = --Vu k+1 . 

T 

• Pressure update: The last sub-step of the algorithm consists of updating the pressure as 
follows: 

(2.5) p k+ i =p k -i + fc+ 5 - x Vu k+ i. 



Remark 2.1. The parameter x ^ m (2.5) is user dependent. By analogy with the projection-based 
pressure correction schemes, we say that the method is in standard form if x — and the method 
is in rotational from if % > 0. 



Remark 2.2. The splitting of the momentum equation in ( 2.2 )-( 2.3 1 is obtained by using the original 
alternating directions ( ADI) scheme of Peaceman and Rachford, see [IB] ■ 



Remark 2.3. The quantity cj)*'~2 can be estimated in many ways. For instance one can take i 
0; this limits the convergence of the scheme to first-order. One can also take </>*'~3 — p*>2 — p 
where is p*'2 is an estimate of p(f). 



A remarkable feature of the algorithm (2.1l to (2.5) is that, although the Dirichlet boundary 
condition on the velocity is not enforced on the entire boundary at the integer time steps, it is 
indeed fulfilled as claimed in the following 



Proposition 2.1. Let {u k } be the velocity sequence from the algorithm (2.1) to (2.5). Then u k \gn = 
for allk = 0,..., K. 
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Proof. It is clear that the boundary condition is satisfied at y = 0, 1. Now taking the difference of 
(2.3) and (2.2), we obtain the following expression for the half-step velocity 

,fc+i i „k 



12.01 „*+J = !L_tIL_^ w(u * + i _„'-,. 



Let us consider x = 1, the other boundary can be treated similarly. The boundary condition at 
x = 1 on the half-step velocity u k+ ^ implies that 

u k+l (l, y) + u k (l : y) = T -d yy {u k+l {l, y) - u k (l, v))- 

Moreover, the boundary conditions on u k+1 and u k at y — 1 imply that this can be re-written into 
the following evolution equation 

u k+l {l, y) - u*(l, y) - T -d yy {u k+1 - u k )(l, y) = -2 M fe (l, y), (u k+1 - u k )(l, -)lv=o,i = 0. 

Since u°(l, y) = and the evolution operator is positive definite, we obtain that u k (l, y) — for all 
k = 0,...,K. □ 

This result turns out to be crucial for the error analysis. 

2.2. Three Space Dimensions. The purpose of this section is to propose a three-dimensional 
version of the above splitting technique. Since the alternating directions method of Peaceman and 
Rachford described in [18] does not extend to three dimensions, we use the alternating directions 
method proposed by Douglas [1] instead to approximate the momentum equation. 

The algorithm is again composed of four steps: pressure predictor, velocity update, penalty step, 
pressure update. 

• Pressure predictor: Denoting by p the pressure field at t = and 0*'~2 an approximation 

of |r9tp(0), the algorithm is initialized by setting = p and <f)~i = </>*'~3 . Then for 
all k > a pressure predictor is computed as follows: 

(2.7) p*' fc+ 5 =p*-S +4> k - 1 2. 

• Velocity update: The velocity field is initialized by setting u° = uo, and for all k > the 
velocity update is computed by solving the following series of one-dimensional problems: 
Find rj k+1 , C fc+1 , and u k+1 such that 



cfc+l _ fc 

S — - A U k + X7n*' k +h = ffc+3 

1 



(2.9) ^ ^ ^d xx (r) k+1 - u k ) = 0, 7 ? fe+1 U =0 ,i = 0, 





T 






-1 _ 






r 






-l 










u k - 


-l _ 





(2.8) 5 ^ _ Au k + Vp*- k +i = f k +i , ^ fe+1 | an = 



l^.iOi — \d yy (( k+1 -u k )=0, C fe+1 | y =o.i = 0, 

(2.11) 5 \d zz (u k+1 - u k ) = 0, U fc+1 U =o ,i=0. 

• Penalty step: The pressure-correction </> fe+ 2 is computed by solving 

(2.12) A4> k+ ^ = --V-u k+1 . 

T 
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• Pressure update: The last sub-step of the algorithm consists of updating the pressure as 
follows: 



(2.13) = j/'-* + K+ i - xV-u* 

Remark 2.4. This method is an extension of the alternating direction method proposed by Dou- 
glas [3]. In order to see this, we add (2.8) and (2.9) to obtain 



V 



k+l 



-\d xx (r] k+1 +u k )-{d yy 



d z 



Adding (|2.8|) ( |2.10[ ) we obtain 

fc+i _ ,,fc 



c 



1 



d xx {r] 



fc+i 



u k ) 



dyy{C k+1 



d zz u k + Vp 



ffc+| 



Finally, adding (2.8)-(2.11) we obtain 

,fc 



u k+l _ u 



1 



1 



- ^d xx (v k+1 + u k ) - ^d yy (c k+1 + u k ) - 



d zz {u 



k+l 



*,fc+i 



pfe+i 



These equations correspond to (3.1a)~(3.1c) of [3], respectively. 



Proposition 2.2. Let {u } be the velocity sequence from the algorithm ( 2.7 ) — ( 2.13 1 . Then u \dQ = 
for allk = 0,...,K. 

(2.14) M fe | ao = 0,Vfc>0. 

Proof. Proceed as in the proof of Proposition |2.1| □ 

2.3. Compatibility Conditions. Note that p := p|* = o is n °t part of the initial data but this 
quantity can be computed by solving 

(2.15) A Po = V(/o + Au ), 9„ Po |r = (/o + Au )-n, 

where we have set fo :— f\t=o- This then requires the initial data to satisfy the following compati- 
bility condition at the boundary (— Auo + Vp — /o)|r = which we assume to hold. This condition 
holds for instance if Uo = and fo = 0, i.e., the fluid is a rest at t = and the source term is zero 
at t = 0. If the above compatibility condition is not satisfied, the error analysis must be adapted 
to account for weighted error estimates by proceeding as in [TH [TS] . 

3. Error Analysis of the Standard Scheme 



The purpose of this section is to study the convergence of the algorithms ( 2.2 )-( 2.5 ) in two space 
dimensions and (2.7)-(2.13) in three space dimensions for x = 0. The main claim of this section is 
that the standard version of our scheme is unconditionally stable and first-order convergent in all 
quantities. 

3.1. Consistency of the Momentum Equation. To evaluate the consistency error on the mo- 
mentum equation, we re- write the momentum equation in a more recognizable Crank-Nicolson form. 
This is done in two space dimensions by adding (2.2) and (2.3) as follows: 



(3.1) 



- d x 



dyyU 



*,fc+i 



= / 



fc+i 



Then using (2.6) we obtain the evolution equation for the integer steps, 
(3.2) u^-u k _ A _ k+ , + Vpk>k+h + r j +l = /fe+i 



DIRECTION SPLITTING 



The same trick can be used in three space dimensions as suggested in [3]. By proceeding as 
above, the intermediate steps, £ fe+1 , r) k+1 , and C fc+1 : can be eliminated, so that the momentum 
equation becomes: 



(3.3) 



u k+l _ u k 



Au k+ ? 



d z 



xyy T vyy ZZ -r v ZZ xx ^xxyyzz) 



:0, 



)5u k+1 + Vp 



Owing to the definition of the operator B (see (1.16l), the momentum equation can be re-written 
as follows independently of the space dimension: 



(3.4) 



u k+l _ u k 



-Bdu 



fe+i 



3.2. Consistency Analysis of the Algorithm. Let u, p be the solution of (1.1 1. We define the 
following velocity and pressure errors: 



(3.5) 



e k+1 := u k+1 - u k+1 . e fc+ 5 



p fc +3 _p*H 



where u fc+1 := u(t k+1 ) and p fc+ 5 := p(t k+ i). 

Next, we obtain equations controlling the errors. Since \ = 0, the pressure update implies that 
the pressure predictor can be written as follows: 



(3.6) e*' fc+ 5 = 2e k -i - e k -$, p*- fe+ s = 2p R 2 _ p - 

that is, the pressure predictor is a second-order extrapolation of the pressure at time level k + | , 
and upon subtracting (3.4) from the momentum equation (1.1 1, we obtain 



(3.7) 



(1 + — B)(e k+1 - e k ) - rAe fe+ 5 + rVe*' fc+ 5 = fR*+* , 



where the residual lZ k+ 2 is defined by 



(3.8) K k+ ^ 



Su k+1 



(u t ) 



fc+i 



V 



,*,fc+i 



~B[5u 



fc+ii 



Finally, using pTH] ) (or pig} ) with x = to eliminate fc +5 from (O) (or pig] )) and using the 
incompressibility constraint, we obtain 



(3.9) 



(5e k +Kq)A = -(e k+ \Vq) + (6p k+ Kq) J 



Vq e Y. 



Note that it is not legitimate to write the equality in strong form, i.e., Ade k+ i is equal to — 'V'e fc+1 + 
A6p k+ i since Sp k+ ^ is not in D(A) (i.e., Sp k+ i does not satisfy the artificial boundary conditions 
associated with D(A)). 

Lemma 3.1. Let u € W 2 -°°{H 2 (n)) n W l '°°{D{B)) and p € T4 72 '°°(i/ 1 (il)). TTien 



(3.10) 



2r(^ fc+ 5, w ) < cr 5 + 7||u||^2, VveL 2 (0). 



Proof. Each of the terms in 1Z k+ i is C(r 2 ), given the smoothness of the exact solution. Note that 

pfc+2 _ p*,fe+5 = <5 2 p fc +2. □ 
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3.3. First Order Estimates on the Velocity. Let us assume that the quantity 0*'2 is estimated 
so that the following holds 



(3.11) 



llP(i) <CT. 



This is the case if (j)*>2 = and if the pressure is smooth enough, say p € C°([0, T], L 2 (fl)). Then, 
the main result of this section is the following first-order convergence statement: 



Theorem 3.1. Assume that the solution (u, p) to (1.1 1 is smooth enough (say u € W 2 '°°(H. 2 (Cl)) n 
W 1 '°°(D(B)) and p € W 2 '°°{Y)). Then, provided that ( |3 11| holds the solution {u T ,p T ) to the 
discrete scheme (2.1) -(2.5) in two space dimensions and (2.7) -(2.13 1 in three space dimensions, 
with x = 0, satisfies the following error estimate 

(3.12) ||e T ||^( L 2) + ||e T ||f2( H i) + r||e T ||^«.(fl) + T\\e T \\ e ^ {A) + v / r||<5e T || £ 2( B ) < cr. 



Proof. Multiply equation (3.7) by 2e k+1 and integrate over ft. Since both the exact velocity and 
the approximate one at integer time steps satisfy the full boundary conditions, we obtain 

(3.13) (l-7)||e fe+1 ||^ + ||6e fe+1 ||£ 2 + J||e fc+1 ||^ + 2T||e fc+ *||f I1 + 2r(Ve*' fc +5, e fc+1 ) + 



\5e 



fc+l||2 

Is 



< e 



fc||2 

II 2 



„fe 1 1 2 



CT 5 . 



Where we have used the identity 2a(a ± b) — a 2 — b 2 + (a ± b) 2 . 
By using 2r 2 e* ,k+ ^ as test function in (3.9) We obtain 



2T 2 (6e k+ i,e* 



2T 2 (5p k+ ^,e^ k+i - 



Clearly, 



-i/x^k+h c *,k+±\ . _ 2 T (e fe+1 , Ve*' fc+ 5 l 

( ( 5 c *+i ) e *' k+ i) A = (Se k+ i ,e k+ i) A - (Se k+ ^ , S 2 e k+ ^) A , 



2 u2 ,f„ u\2 



so that using again the identity 2a (a — b) = a — b + (a — b) we obtain 



(3.14) r 2 ||e 



,:+l lli-l|e fc ^ni 



We* - » 



.2|U2 fc+i||2 



2r(e 



fc+i 



Vr 



*,fc+i> 



2r 2 (,5p 



A- 



To obtain a control on ||£ 2 e fc+ 2 H^, we apply the time increment operator 6 to (3.9) (assuming that 
k > 2) and we use the test function rS 2 e k+ i: 



So that 
(3.15) 



\\S 2 e k+ i ||i = (Je fc+1 , V5 2 e fc +5) + r(<5 2 p fe +5 , S 2 e k +i) A 

< ||5e fc+1 || L 2||V<S 2 e fc+ 3|| L2 + r^p^ |U||<5 2 e fc +i 



< (||<5e fc+1 || L 2+r||5 2 p fc+ S|U)||^e 



.2 fc+i | 



-^ 2 || ( 5 2 e fc +i H2, < H^e^^ 1 !!^, + ^ 2 ||^ 2 p fc +* ||2 ^ 2^||5 2 p fc +^ lUH^e^^ 1 !^.. 
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Adding (|3l3|), (|3T4|) and (|3J5]) we obtain 

2r||e 



(l-J)ll. 



fc+l||2 
lL 2 



2 l|e 



fc+l 1 1 2 



0|e fe+1 ||| + WSe^Wl] <cr 



H 1 

fc+l |2 



fc+i||2 
H 1 



|e fe ||^ 



«H-i||a 



\5e 



Hi- ill 2 



H 1 



2 



+ 2r 2 (<Jp fe+ 5, e*< k+l i) A + T 2 \\5 2 p k+l i \\\ + 2r\\5 2 p k+ i 

Let us examine the last three terms in detail: 
• T 2 ||(5 2 p fe+ 2 11^. Given the smoothness of p this term is C(r 5 ). 
*2„fc+5ll 1 ||(5e fc+1 ||L2. We estimate it as follows: 

„fc 



l|e fc | 



fc||2 
B 



\8e 



fc+i i 



|L- 



• 2T\\S 2 p k ■ 2 \\ A \ 

2r||5 2 p' £ +*|U|| ( 5e' £+1 || L2 <cr 3 (|| e ' £ + 1 || L2 



3 fc+l||2 

Hl 2 



? k \\h- 



• 2t 2 (ASp k+ i , e*' k+ i). Given the smoothness of p 



2T 2 (Sp k+ 3, e*> k+ i) A = 2r 2 (<5p fc+ 3 , e k -^) A + 2r 2 {5p k+h -, 6e k ~^) 

■cT 3 \\6e k -^\ 



<cT 3 \\e k - l *\\ A 



< CT 6 



r 3 ||e fc -^ 



2 

\A 



r 2 \\Se 



A 

fc-f 



|2 
\A- 



Note that this term is the only one in the entire error analysis that spoils the game. This consistency 
term does not allow us to obtain directly an error estimate of order larger than O(t). 
We have finally proved that the following holds for all k > 2: 



(1 - -)||e fc + 1 ||^ 2 + -| 



= fc+l||2 

IIh 1 



\Se k+1 \\ 2 B ] <cr 3 + (l 



2r\\e k ^\\ 2 H 

T 



2||,fc+i 1 1 2 

u 



fc||2 



t ||e 

T 



[II 



fc+l||2 

lis 



„fc||2 



IL 2 ^-^ IIhi+T" + 

Upon observing that the initialization process (p~^ = p ) implies 

r 2 \\S 2 el\\ 2 A <(l + r)\\e 2 \\ 2 +r^ 



fe ^lli 



„fc||2 



we infer that the above inequality holds also for k = 1. As a consequence of (3.11), we also deduce 
that 



\e'\\h 



Hi 



l H1 



Vile 1 II 2 <cr\ 



By summing the above relation from k = 1 to K and by applying the discrete Gronwall lemma 
allows us to conclude. □ 

The ability of 8u k+1 /r to approximate u t is made explicit in the following: 



Lemma 3.2. Let the solution (u,p) to (11} be smooth enough u e W 3 '°°(H 2 (fi)) n W 2 '°°(D(B)) 
and p € W 3 '°°(Y)). Then the solution (u T ,p T ) to the discrete scheme (2.1 )-( [275| ) in two space di- 
mensions and (2.7)- |2T3 | three space dimensions, with \ — 0, satisfies the following error estimate 

(3.16) ||£e r ||^(L 2 ) + ||<5e T ||^2( H i) + T\\5e T \\ e ^^ A) + r||(5e T || f oo (B) < cr 2 . 

Proof. Apply the arguments in the proof of Theorem |3.1| to the time increments. □ 
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3.4. Error Estimates on the Pressure. It is known that for the incremental projection scheme 
in standard form it is possible to prove that the error on the pressure in the £ 2 (L 2 )-norm is C(t) 
(cf. [3 [H[22]). The purpose of this paragraph is to show that, although on a weaker norm, a similar 
result holds for the proposed algorithm. Let us define the norm 

<yq,v) 



(3.17) 



MIa 



sup 
o^cez \\v\\z 



Theorem 3.2. Assume that the hypotheses of Lemma 3.2 hold, then 
(3.18) IM^A) < cr. 



Proof. Using the error equation (3.7) we obtain 

1 



|e*' fc +!|| A 



SUP T-T- 

O^vEZ \\V\\Z 



^— , v) + (Ve fc+ i , Vw) + - (Se k+1 . 



< 



\\Se k+1 \ 



L 2 



|e fe+ ^|| H1 



Tll^ fe+1 || B +cr 2 



< cr - 



|e fe ^|! Hl 



where the last estimate holds in view of Lemma 3.2 Take the square of this inequality, multiply it 
by t and sum over k. The result follows by using the conclusion of Theorem |3.1| □ 



Remark 3.1. It seems that it may be possible to obtain a first-order error estimate on the pressure 
in the £ 2 (i 2 )-norm in the fully discrete case under the additional (somewhat restrictive) condition 



(3.19) 



r < 




This is a CFL condition in two space dimensions. The reasoning behind this conjecture is the 
following. Assume that the velocity is approximated using a finite-dimensional space and that 
the norm in B is appropriately approximated, say || • \\b h - In view of (1.181 it is reasonable to 
expect that the following inverse inequalities hold: 



IMk < 



-i 



ch "IIUIIhi 



c/l -1 (l + T2h" 



u|| H i m 



Then, assuming that the pressure is approximated using a space C H^ (Q) so that the pair 
(X/^Mh) satisfies the so-called LBB condition, [S1|B], we obtain 



\ L i < sup 



< sup — n — 

k + 1 



(6e k+1 /r,v) + (Ve fe+ 5, Vv) + T -{5e k+ \v) Bh + (K k+ Kv) 



I^IIh* 



The two-dimensional inverse inequality implies 
\\e*> k+ i\\ L 2 <cr 2 + \\T- 1 Se k+1 \\ Ij , 



T 

— sup 
4 o/i>ex h 



\Se k+1 \ 



| h i +cr 2 ^- 2 ||r- 1 ( 5e fe+1 ||Hi 
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whereas the three-dimensional inverse inequality implies 

||e*' fe +i || L2 < cr 2 + \\r- l 5e k+1 \\v + ||e fc+2 1 H1 + c(t 2 ^ 2 + T*h- i )\\T- 1 5e k+1 \\ w i. 

Take the square of this inequality, multiply it by r and sum over k, then the estimates of Lemma |3.2| 
together with condition (3.191 yield the desired estimate, ||e T |U 2 (i, 2 ) < cr. 

3.5. Second-Order Estimates on the Velocity. Despite the fact that numerical experiments 
suggest that the standard form of the above algorithm is close to second-order on the velocity in the 
L 2 -norm, (see Section [6|, a proof of such statement eludes us at the moment. We briefly elaborate 
in this section on the difficulties that arise when trying to establish a second-order error estimate. 

The argument one usually invokes to prove a second-order error estimate consists of multiplying 
the error equation by Se k+ i, where S is the right-inverse Stokes operator (see (1.22)). Following 
this reasoning, and using property (1.24), we obtain that the following holds 



(3.20) i (|e fc+1 | 2 - |e fc | 2 ) + ^ \\e k+ i || 2 2 + ^(<5e fc+1 , Se k +^) B < r(n k +K Se k+1 ) 

+ cT\\e k+1 i -P a e k+ i\\l 2 , 

Provided the exact solution is smooth enough, we can estimate the residual term in a way similar 
to Lemma 13. 11 



T(K k+ KSe k+ ^} < ct 5 



2 

L 2 ■ 



Using the estimates of Lemma 3.2 we can control the i?-norm as follows: 



— (Se k+i ,Se 



< -^\\Se k 



fc+ii 



\Se k+ h\ B < cr 3 ||5e fe+ !| 



In two space dimensions the H 2 -regularity of S implies ||Se fc+ 2 || B < c||e fe+2 || L2 so that 



— (5e k+i ,Se 



fc+i||2 

II 2 - 



Note that the above reasoning does not apply in three space dimensions. In conclusion, in two 
space dimensions (3.20) becomes 

\ a k+l\2 



e k \l + T\\e k+ Hl. < ct(t* + inf ||g fc+ * - v\\i 2 ) = ct(t* + \\e k+ i P u e k+ i ||£ a ). 



which in turn yields 
(3.21) 



||e T ||| 2(L2) < c(r 4 + \\e T - P H e T || 2 2(L2) ). 



This inequality shows that the estimate on ||e r ||^ 2 (L 2 ) is controlled by \\e T — Pne T ||£ 2 (l 2 )- Let us 
now try to bound |je r — Pnc T \\e 2 (l 2 ) uniformly. 

By definition, there is /i fe+ s 6 H^ =0 (tt) so that e fe+ 2 — Pne k+ ^ = W/J, k+ ^. In other words fi k+ i 
solves — A/i fc+ 3 = — V-e fe+ 5 and d n fj, k+ ^\gn = 0. Then the penalty equation (3.9) together with 
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the assumed smoothness of the pressure and the estimates of Lemma |3.2| imply that 



Pne k+l i\\h 



\V^\\h 

<£(l|fc* + *IU- 

<cr 2 \\^\\ A . 



fe+i 



Se k 2 - 5p k 



This finally gives the estimate 



P H e fc+ 5||L2 <cr 2 - 



Im* + *IU 



A induces a norm equivalent to H 1 . This is unfortunately 



which can be controlled uniformly if 
not true with the operators A defined in and ( 1.12 1. 

In conclusion, the reasoning carried out above seems to indicate that the right-inverse Stokes 
operator S is not the correct operator that should be used for the duality argument. The operator 
that should be used instead still eludes us at the moment. 



4. Error Analysis of the Rotational Scheme 



The purpose of this section is to analyze the algorithms ( 2.1 H 2.5 1 and ( 2.7 )-( 2.13) for \ ^ 
and to show that, as it is the case for the classical rotational pressure-correction schemes (cf. [13]). 
these algorithms provide a better order of convergence than the standard form. 



4.1. Consistency Analysis. Let u, p be the solution of (1.1). We define the following velocity 
and pressure errors: 



(4.1) e k+1 := u k+1 e"" =» := p"' => - p" 

where u fc+1 := u(t k+1 ) and p fc +i := p(t k+ i). The error on the pressure correction is measured by 
introducing the following quantity: 



,fc+i 



(4.2) 



Using the above notation we infer 



P 



k 3) = p fc+ 5 - (p fe ~3 



= P" ' 2 - (P 

= 5 2 p k+ ^ + (e fe_ 5 + 5e k ^^ + xV-e fc ~5") 
Then momentum equation is rewritten as follows: 



X V-u k -i) 



(4.3) 



k+l 



e fc ) - rAe* 



fc+i 



rV(e 



fe-7 



) = TK 



where the residual TZ k+ = is defined by 



(4.4) 



Su k+1 



"(ut) 



A 



- v 



-B[6u k 



fc+n 
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The result of Lemma 3.1 holds again, lZ k+ 2 = 0(t 2 ), provided the exact solution is smooth enough. 
The equation that controls the pressure correction is re-written as follows: 

(4.5) (V k+ Kq)A = -(e k+ \Wq) + (Sp k +^q) A Vq E Y. 

T 

4.2. A Priori Estimate on the Divergence of the Velocity. Let us assume the that quantity 
is estimated so that the following holds 

(4.6) ||P(5)-J»*' 4 IU» <cr. 

This is the case if <fi*'i — and if the pressure is smooth enough, say p € C°([0, T], L 2 (f2)). The 
main result of this section is the following 



Theorem 4.1. Assume that the solution (u,p) to (1.1 1 is smooth enough, (say u € W 2 '°°(H 2 (fi)) fl 
W 1 '°°(D(B)) and p G W 2 '°°(Y)). Then, provided (4.6) holds, the solution (u T ,p T ) to the discrete 
scheme (2.1 ) — ( 2.5 ) in two space dimensions and (2.7)-(2.13) in three space dimensions, with < 
X < I, satisfies the following error estimate 

(4.7) ||<W|| 2 oo (L3) + t||Vx<M^ (l2 



|Vx<5e r " 2 



£ 2 (L 2 ) 



rllVeJ 



< CT 4 . 



Proof. Following |13j . we derive an improved estimate on the divergence of the velocity. This is 
done by working with the time increments of (4.3)-(4.5). 

Apply the time increment operator S to the momentum equation (4.3) and test against 28e k+1 
to obtain 



(4.8) 



1 - 



H^+i 2 , + \\S 2 e^\\l + T -\\6e k +^ + 2r\\6e k ^\\ 2 H1 + -\\5e k +^ 2 



B 



2r < V [Se 



SV k 



, Se 



k+i 



< CT* 



\\Se k \ 



T 



2ll<fe*llHi + \\ 



6e k \\% 



where we used the fact that the residual is C(r 2 ). Note that we could decrease the consistency error 
to C(t 3 ) by assuming more regularity on u and p, but it would not improve the overall accuracy of 
the method since the splitting error will turn out to be 0(t 2 ) (see below). 

Apply the time increment operator S to (|4.5|) and use 2r 2 V k+ 2 as a test function. We obtain 



r 2 \\V 



k+i it 2 



T 2 \\5V k+ * 



2 

A 



T 2 \\V k 



-2T(W-Se k+1 ,V k+ ^) +2T 2 (S 2 p k+ i,V k ■ -) A 
= -2T{V5e k+1 , X Ve k+ i + Se k+ ^) + 2r 2 (S 2 p k+ ^ , V k+ ^) 



where we used the identity 2a(a ±6) = a 2 + (a ± b) 2 — b 2 . This gives 

(4.9) T 2 \\V k+ ^\\+T 2 \\SV k+ ^\ 2 A + X r\\^e k+1 \\h = r 2 \\V k -^\\ 2 A 

+ xr||Ve fe || 2 2 + 2 T (6e k+1 , VSe k+ ^] 



2T 2 (S 2 p k+ i,V 



Apply again the time increment operator S to ( |4.5| ) and test with — 2T 2 5 2 T >k+ ^ . Using, again 
the identity 2a(a — b) = a 2 + (a — b) 2 — b 2 , we obtain 



(4.10) -T 2 \\\5T k+1 2\\ 2 A + \\5 2 T k+ ^\\ 2 A -\\5V k -^\\ 2 A \ = -2r(VS 2 r k+ K 5e 



k+l\ 



2r 2 (S 2 p k+ K5 2 V k+i - 



Observe that (4.9 ) + (4.10 1 amounts to testing the time increment of (4.5 ) with 2t (V 2 -\-SV 2 ' 
We have split the two steps to make the argument clearer. 
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By summing (4.8), (4.9) and (4.10) we deduce that 

(4.11) (l - I) \\6e k ^\\h + \\\6e^\\^+2r\\5e k ^\\\ l , + + xr||Ve fc+1 || 2 L2 

+ r 2 \\V k+ i\\ 2 A +r 2 \\6V k -if A + \\6 2 e k+1 \\l* ~ T 2 \\S 2 V k+ *\\ 2 A - 2 X r(VSe k +i ,VSe k+1 ) < ct 5 
+ \\Se k \\li + ^11^11^ + ^||5e fe |||+xr||Ve' £ ||| 2 +r 2 \\V k -i \\\ + 2t 2 {5 2 ^ ,V k ^- + 5V k -^) a- 
where we used the following identities 

e fe-5 + 3V k -^ - 5e k+ i + S 2 V k+ % = X VSe k+ K 

■pk+i _ fi2rpk+± _ -pk-1 _|_ ftpk-\ 

Given the smoothness of p, the following holds: 

3 

2T 2 (S 2 p k+ ^V k -^ +SV k -^) A < CT 5 + y||P fe "5 \\\ +^11^-2 

Observe that it is here that the irreducible splitting error comes into full light. Although the 
consistency of the time increment of the momentum equation is 0(t 3 ) (provided enough regularity 
is assumed on u and p), the above inequality shows that splitting error of the method is 0(t 2 ). 
Then (4.11) becomes 

(4.12) (l-I) Pe^l^ + T -\\6e k +Yni ~ x\\\^e k+l \\ 2 L , +2r||^^||^ - 2r X || || 2 2 

+ l\\Se k+1 \\ 2 B + XT||V-e fe+1 ||| 2 + T 2 \\P k +i \\ 2 A + ||<5 2 e fe+1 1| 2 2 - r 2 \\S 2 V k+ ^ \\ 2 A < cr 5 

+ \\Se k \\h + p^llH! - lx\\VSe k \\h + Jl^lll + XT|| Ve & ]|| 2 + r 2 (l + T -)\\V k ^ f A . 
where we used 

T „„ „ u i_ 1 mo T 



~2 X (VSe k+ Kv5e k+1 ) = - X II V<5e fc+1 1| 2 2 - -|| V<5 e fc ||| 2 + 2r\\ We fc +^ \\ 2 L , 
Then the identity §TJ\ gives 

(4.13) (l - J) ||^ fc+1 || 2 2 + pVx<5e fc+1 || 2 2 + (1 - x )l\\V.Se k+1 \\ 2 L2 + 2r\\VxSe k +i 
+2^(l- X )||V-5e fe +* ||i 2 + J||^e fc + 1 ||2 ?+x ^|| V . e fc + 1 ||i 2 +^ 2 ||7^ fc +i||^+||^ 2 e fc + 1 ||2 2 _^2|| 5 2^A ; +i ( j ^ ^, c ^ 5 

+ ll<Se fc || 2 2 + pVx5e fc || 2 2 + (1 - x)pWe fe ||! 2 + ^||,5e fc ||| + xr||V- e fc || 2 2 + r 2 (l + T -)\\V k ^ \\\. 

To conclude we are going to observe that the quantity ||<5 2 e fc+1 ||l 2 — T 2 \\5 2 V k+ ^ \\ 2 A is non negative 
up to some consistency error. To see this, let us apply the time increment operator 5 2 to (4.5) and 
test the equation with TS 2, P k+ ^ . After using the Cauchy-Schwarz inequality and the inequality 
(1.8), we obtain 

r\\S 2 V k+ ^\\ A <\\6 2 e k+1 \\ L2 +r\\S 3 p k ^\\ A , 
which, given the smoothness assumption on p, then implies 

T 2 \\8 2 V k+ H 2 A < ||^ e fc + 1 || 2 2 +r 2 || ( 5 3 p' £ +5|| 2 4 + 2r|| ( 5 2 e fe + 1 || L2 || ( 5 3 p fe+ 5|| A , 



<cr 5 + \\6 2 e k ^\\h + \\\5e k +'\\l 2 + T -\\&e k \\l 2 . 
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Note again that the consistency error could be decreased to C(r 3 ) by assuming p <E W 3 '°°(Y), but 
this would be useless since the splitting error of the method has been shown to be 0{t 2 ) above. 
By adding this last inequality to (4.13) we finally obtain that the following holds for all k > 2. 



1 - 



\Se 



fc+l||2 



|Vx5 e fc + 1 ||2 2 + (l- x )I||V5e fc+1 ||| 2 +2r||Vx5 e - fc +5||2 2+ 2r(l~ X )||V-<5e 



2 

L 2 



\Se 



fc+l||2 



L 2 
T 



b + XT||V-e*+i£ 2 + T 2 \\V k+l2 \\\ < cr 5 + (1 + ^)\\Se k \\h + ^||Vx<5e fc ||^ 



(1 - x)-||We fc ||i 2 + -\\5e k r B + xr||V-e fc ||i 2 + r\l + 



li. 



The following estimate also holds as a consequence of the initialization hypothesis (4.6 1: 



\5e 2 \\h + r\\Vx8e 2 \\l 2 + r||V-,5e 2 ||| 2 + T\\5e 2 \\ 2 B + X^I|V-e 2 ||| 2 + T 2 \\pi | 



< r 4 . 



By summing the above inequalities from k = 2 to K and by applying the discrete Gronwall 
lemma we finally obtain the following error bound: 



||<5e r ||^L (L2) +T||Vx,5e T | 
This completes the proof. 



^(L 2 



r|| Vx5e 1 -||| 2 ( L 2) 



XT||V-e T |||ao (i2) 



\Wr\ 



e°°(A) 



□ 



4.3. Error Estimates. Having obtained the estimate of Theorem 4.1 we can now show that the 
rotational version of the algorithm provides a better order of convergence for the velocity in the 
£ 2 (L 2 )-norm, at least in two space dimensions. To this end, let us denote by ip T the sequence whose 
generic term is ^ fc+2 := ^(ifj k+1 + ip k ) 

Theorem 4.2 (^ 2 (L 2 ) Velocity Estimates). Assume that the space dimension is two. Under the 
assumptions of Theorem 4-1 the solution (u T ,p T ) of the scheme (2.1 )-(2.5) in two space dimensions 
satisfies 

|| U T - W T ||^(L 2 ) < CT?. 

Proof. The proof proceeds by a duality argument using the right-inverse Stokes operator 5". By 
proceeding as in Section 3.5 we obtain (see (3.21)): 

||e T |l^(L2) < c(r 4 + ||e T - fke r ||| 2(L2) ). 



The estimate (4.7) then immediately implies 



IM? 2 (l 2 ) < c(r 4 + ||Ve r || 2 2(L2) ) < cr 3 . 



Now we observe that 



= fe+l||2 



\Se k+ 



1 II 2 
L 2 



< 



!^|| 2 2 



ll^ +1 HL, 



which, along with (4.7), implies 



ll e "H|^(L 2 ) 

which completes the argument. 



< 



lerll^fL 2 ) +c||<5e fc+1 || 2 00(L2) <cr 3 , 



□ 



Let us now show convergence of the velocity in the € 2 (H 1 )-norm without any restriction on the 
space dimension. 



Theorem 4.3 (^(H 1 ) Velocity Estimates). Under the assumptions of Theorem J^.l the solution 
(u T ,p T ) of the scheme (2.1)-(2.5) in two space dimensions and (2.7)-(2.13) in three space dimen- 
sions satisfies 



= (Hi) 



< CT. 
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Proof. Observe first that the following holds for all k £ {0, 

k 

||Vxe fc +5|| L2 < ^||Vx«5e l+ 5 
»=i 

which implies 



,K}: 

|Vxe5|| L 2, 



||Vxe fe+ 5|| L 2 < ct~ 1 \\ Vx<5e T ||f2 (L 2 ) + ct < c(r V 2 + r) < cr, 
and, owing to (1.7 1, this concludes the proof since we have already established that || V-e k+ i \\ L 2 < 



a 



Remark 4.1 (Pressure Error Estimates). The same methods and ideas used in Section 3.4 can be 
invoked to show that the pressure satisfies the following estimate 



U 2 (A) 



< CT. 



We omit the details for the sake of brevity, 



Remark 4.2. Whether Theorem 4.2 holds in three space dimensions is not clear. The main obstacle 
in the way is the splitting error induced by the splitting of the momentum equation. Based on our 
num erical experiments, we conjecture that both the error estimates in Theorem |4.2| and Theorem 



4.3 



can be improved by a T2 factor irrespective of the space dimension. 



5. Other Time Marching Techniques 



As mentioned in Remark 2.4 the velocity update (2.8 1— (2.111 is a sequence of three approxi 



mations of the momentum equation where each approximation consists of evaluating the second 
derivative in one of the spatial directions implicitly with the Crank-Nicolson scheme whereas in 
the other directions it either employs the solution from the previous time level, if no implicit 
approximation is yet computed in the given direction, or uses the already computed implicit ap- 
proximations. This observation leads us to propose the following split version of the second-order 
backward difference scheme (BDF2) to approximate the momentum equation: 
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We now write the full BDF2 algorithm in a form similar to (2.7)-(2.13). To simplify the presen- 



tation, let us assume that proper approximations of the velocity and the pressure time derivative 
are available at t = —t and t = 0. If these quantities are not available, we start the scheme with a 
lower-order approximation at the first time step in order to compute those approximations. 

• Pressure predictor: Denoting by p the pressure field at t = 0, by 0*'° an approximation of 
_1 an approximation of rd t p(— r) the algorithm is initialized by setting 
and (j)^ 1 — (f)*^ 1 . Then for all k > a pressure predictor is computed 



rd t p(0), and by </>*' 
P " = Po , 0° = P>°, 
as follows: 



(5.1) 



p *,k+i =p k 



k-l 
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Velocity update: The velocity update is computed by solving the following series of one- 
dimensional problems: Find £ fe+1 , r/ k+1 , £ fc+1 , and u k+1 such that 



(5.2) 
(5.3) 
(5.4) 
(5.5) 

(5.6) 



4u k + u k 
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3(C fc+1 - 




2r 




3(u k+1 - 
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~d zz (u k 

Penalty step: The pressure-correction 



-d yy (( k+1 -u k ) = o, c k+1 \y=o,i = o, 

k+1 -u k )=Q, u k+1 \ z=Q:1 =0. 

fe+1 is computed by solving 
3 



2r 



fe+i 



• Pressure update: The pressure is updated as follows: 
(5.7) P k+1 = P k + <f> k+1 - xV-u fc+1 . 

Note that this scheme is formally second-order consistent because eliminating the intermediate 
velocities results in a second-order perturbation of the classical pressure-correction BDF2 scheme. 
Numerical experiments show that this algorithm is indeed unconditionally stable when tested on 



the unsteady Stokes problem and its rate of convergence is similar to that of (2.7 1— ( 2.13 ) . 

6. Numerical Experiments 

We report in this sections numerical tests aiming at evaluating the performance of the algo- 
rithms (2.1 )— (2.5) in two space dimensions and (2.7)-(2.13) in three space dimensions. The space 
approximation is done using the MAC scheme. 



6.1. Accuracy Tests. The standard and the rotational versions of the scheme (2.1 1— ( 2.5 ) have 
been tested numerically on a two dimensional analytic solution and the results have been reported 
in [9]. The rate of convergence with respect to r for the velocity in the L 2 -norm for both versions 
of the method is about 1.8 or higher, whereas for the pressure in the L 2 -norm it is about 1.85 for 
the rotational version and about 1.5 for the standard version. 

We now investigate the convergence rates in three space dimensions in fl = (0, l) 3 using the 
following solution of the unsteady Stokes problem (with the appropriate source term) : 

U! = (sin x cos y sin z — sin x sin y cos z) sini 

U2 — (sin x sin y cos z — cos x sin y sin z) sin t 

U3 = (cos x sin y sin z — sin x cos y sin z) sin t 

p = cos(a; + y + z + t). 

We display in the left panel of Figure [l] the L 2 -norm of the error on the velocity at T = 2 versus 
the time step r for the rotational scheme with x = 1. The L 2 -norm of the error on the pressure is 
displayed in the right panel of the figure. The convergence rate on the velocity varies between 1.6 
and 1.8 while the convergence rate on the pressure is comprised between 1.5 and 1.7. From tests 
not reported here, we have observed that the standard version of the scheme has a convergence 
rate between 1.6 and 1.7 for the velocity and a convergence rate between 1.25 and 1.4 for the 
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pressure. These results suggest that the actual convergence rates of both schemes are higher than 
those theoretically estimated above. However, at the present it is unclear how to improve these 
estimates. 




Figure 1. Rotational form (x = 1). Left: L 2 -norm of the error on the velocity 
(dashed line) at T = 2 on a uniform grid, 100x100 x 100; Right: L 2 -norm of the 
error on the pressure (dashed line). 



6.2. Splitting vs. Projection. To further illustrate the convergence properties of the present 
schemes we now compare it with its unsplit pressure-correction counterpart, i.e., the momentum 
equation is unsplit and the pressure correction is computed by solving the Poisson problem (A = 
— Ajv). The comparison is done in two space dimensions in Q = (0, l) 2 on the following analytical 
solution: 

(6.1) u = (sinxsin(y + t), cosxcos(y + t)), p = cos xsin(y + t). 

We show in Figure [2] the error on the velocity and the pressure as functions r for the unsplit 
second-order projection and the corresponding results using the present direction splitting schemes. 
Clearly, both the standard and the rotational versions of the direction splitting schemes produce 
results that are very similar to those produced by their unsplit counterpart. The largest differences 
are observed on the velocity for the standard version of the schemes. But, even in this case, the 
direction splitting produces errors which are only between 1.2 and 2 times larger than the errors 
produced of the classical standard scheme. The computational complexity of the present schemes, 
however, is significantly lower. 

6.3. Lid Driven Cavity. We compare in this section the performance of the direction splitting 
algorithm with its unsplit pressure-correction counterpart on the so-called lid driven cavity. The 
computational domain is Q — (0, l) 2 . The boundary conditions are u\ x= o,i,y=o = 0; u \y=i = 1 
and v\gn = 0. The computation is done at Reynolds number R e = 100 on a MAC grid composed 
of 40 x 40 nodes and with time step r = 0.01. The advection term is computed by means of the 
explicit second-order Adams-Bashforth approximation. The comparison between the two codes is 
done at t — 1 and t = 10. 
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Figure 2. L 2 -norm of the error on the velocity (top) and pressure (bottom) at 
T = 2 on a uniform 40x40 grid; Left: unsplit projection scheme in standard form 

(dash-dotted line) . Right: unsplit 



(dashed line) and scheme (2.2)-(2.5) with \ = 
projection scheme in rotational form with \ 
([275]) with x = 1 (dash-dotted line). 



1 (dashed line) and scheme ( 2.2 ) - 



We show in figure [3] the horizontal and vertical profiles of the velocity alongside the verti- 
cal/horizontal lines through the center of the cavity. The results of the two schemes (unsplit and 
split) are very close to each other; detailed examination (not reported here) shows that the two sets 
of results differ in the fourth decimal digit. For comparison, we also display with o symbols the 
result of the split scheme on a MAC grid of 200 x 200 nodes and with time step r = 0.0025. The 
three sets of results arc visually indistinguishable. 

6.4. Backward Facing Step. Finally, the new direction splitting method is validated on the two- 
dimensional flow over a backward-facing step. Extensive experimental and computational data on 
this flow is available in pQ and |15j . Here we compute the solution to this problem in a rectangular 
cavity of size 1x16 with a uniform grid of size h — 0.005 and a time step r = 0.001. We prescribe the 
fully developed parabolic profile with maximum velocity | at the upper half of the inflow side and we 
prescribe the no-slip condition at the lower half. At the outlet we impose zero-Neumann conditions 
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Figure 3. The horizontal (resp. vertical) profiles of the velocity alongside the 
vertical (resp. horizontal) lines through the center of the cavity at R e = 100 and 
at t = 1 (left panel) and t = 10 (right panel) on 40x40 MAC grid, r = 0.01. 
Unsplit projection scheme (solid line); direction splitting scheme (□, □ symbols); 
direction splitting on 200x200 MAC grid with r = 0.0025 (o, symbols). 
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on the velocity and the zero Dirichlet condition on the pressure. One important characteristic of 
the flow is the length of the recirculation zone behind the step, say r. We report in Table 6.4 the 
results of the present computations at Reynolds numbers (based on the channel height) R e = 100, 
200 and 400 and we compare these results with those from [15] . The present scheme yields results 
which arc in a very good agreement with the existing data. 



Re 


r/ 


s 




Current result 


Result in 15 


100 


3.22 


3.2 


200 


5.33 


5.3 


400 


8.6 


8.6 



Table 1. Flow over a backward- facing step. Re-attachment length r divided by 
the step height s as a function of the Reynolds number R e for the present compu- 
tations and for the computations of Kim and Moin [15] . 



6.5. Parallel Implementation. We have implemented a parallel version of the algorithm (2.7)- 
( 2.13 ) with the MAC stencil using central differences for the first- and second-order derivatives. The 
algorithm has been implemented in parallel on a Cartesian domain decomposition using MPI. All 
the one-dimensional linear systems are solved in parallel using direct solves of the Schur complement 
induced by the domain decomposition. We have verified that the weak scalability of the code is 
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quasi-perfect up to the maximum number of processors that were available to us without special 
request for allocation, i.e., 1024 processors. 

Extensive numerical tests have shown that the algorithm is stable under CFL condition in the 
Navier-Stokes regime. We have computed a highly accurate benchmark solution for the start-up 
flow in a three-dimensional impulsively started lid-driven cavity of aspect ratio 1x1x2 at Reynolds 
numbers 1000 and 5000. Successive refinements have shown that the velocity field is four digit 
accurate at R e — 5000 for dimensionless times t = 4, 8 and 12. The computations have been 
done in parallel (up to 1024 processors) on adapted grids of up to 2 billion nodes in three space 
dimensions. All these numerical experiments are reported in a forthcoming paper. 

6.6. Further Developments. We believe that the algorithm presented in the present paper has 
a lot of potential for further developments; we are thinking in particular of academic problems 
that can be solved in simple geometries with regular grids e.g. simulation of turbulent flows in the 
atmosphere and in the ocean, simulation of multiphase flows, stratified flows, variable density flows, 
combustion, solution of subgrid problems as part of an homogenization procedure, etc. 

As described in the present paper, the algorithm is suitable only for simply-shaped domains. 
However, there are possibilities to impose boundary conditions either via penalty methods, or 
fictitious domain techniques, or via directional adjustment of the grid at the boundary. The authors 
have implemented the directional adjustment procedure and have observed that the resulting scheme 
is unconditionally stable and convergent for the time-dependent Stokes problem. These results will 
be reported elsewhere. 

References 

[1] B.F. Armally, F. Durst, J.C.F. Pereira, and B. Schonung. Experimental and theoretical investigation of 

backward-facing step flow. J. Fluid Mech., 127:473-496, 1983. 
[2] A.J. Chorin. Numerical solution of the Navier-Stokes equations. Math. Comp., 22:745-762, 1968. 
[3] Monique Dauge. Stationary Stokes and Navier-Stokes systems on two- or three-dimensional domains with cor- 
ners. I. Linearized equations. SIAM J. Math. Anal, 20(l):74-97, 1989. 
[4] Jim Douglas, Jr. Alternating direction methods for three space variables. Numer. Math., 4:41-63, 1962. 
[5] A. Ern and J.-L. Guermond. Theory and practice of finite elements, volume 159 of Applied Mathematical 

Sciences. Springer- Verlag, New York, 2004. 
[6] V. Girault and P. -A. Raviart. Finite Element Methods for Navier-Stokes Equations. Theory and Algorithms. 

Springer Series in Computational Mathematics. Springer- Verlag, Berlin, Germany, 1986. 
[7] J.-L. Guermond. Un rcsultat dc convergence d'ordrc deux en temps pour l'approximation des equations de 

Navier-Stokes par une technique de projection incrementale. M2AN Math. Model. Numer. Anal, 33(1):169- 

189, 1999. Also in C. R. Acad. Sci. Paris, Serie I, 325:1329-1332, 1997. 
[8] J. L. Guermond, P. Minev, and Jie Shen. An overview of projection methods for incompressible flows. Comput. 

Methods Appl. Mech. Engrg., 195(44-47) :6011-6045, 2006. 
[9] J.-L. Guermond and P. D. Minev. A new class of fractional step techniques for the incompressible Navier-Stokes 

equations using direction splitting. Comptes Rendus Mathematique, 348(9-10):581 - 585, 2010. 
[10] J.-L. Guermond and L. Quartapelle. Calculation of incompressible viscous flows by an unconditionally stable 

projection FEM. J. Comput. Phys., 132(l):12-33, 1997. 
[11] J.-L. Guermond and A. Salgado. Error analysis of a fractional time-stepping technique for incompressible flows 

with variable density, submitted to SIAM J. Numer. Anal, 2009. 
[12] J.-L. Guermond and Abncr Salgado. A splitting method for incompressible flows with variable density based on 

a pressure poisson equation. Journal of Computational Physics, 228(8) :2834 - 2846, 2009. 
[13] J. L. Guermond and Jie Shen. On the error estimates for the rotational pressure-correction projection methods. 

Math. Comp., 73(248):1719-1737 (electronic), 2004. 
[14] John G. Heywood and Rolf Rannacher. Finite element approximation of the nonstationary Navier-Stokes prob- 
lem. I. Regularity of solutions and second-order error estimates for spatial discretization. SIAM J. Numer. Anal., 

19(2):275-311, 1982. 



24 J.-L. GUERMOND, P.D. MINEV, AND A.J. SALGADO 



J. Kim and P. Moin. Application of a fractional-step method to incompressible Navicr-Stokcs equations. .1. 
Comput. Phys., 59(2):308-323, 1985. 

T. Lu, P. Ncittaanmaki, and X.-C. Tai. A parallel splitting up method and its application to Navier-Stokcs 
equations. Appl. Math. Lett., 4(2):25-29, 1991. 

T. Lu, P. Neittaanmiiki, and X.-C. Tai. A parallel splitting-up method for partial differential equations and its 
applications to Navier-Stokes equations. RAIRO Model. Math. Anal. Numer., 26(6):673-708, 1992. 
D. W. Pcaccman and H. H. Rachford, Jr. The numerical solution of parabolic and elliptic differential equations. 
J. Soc. Indust. Appl. Math., 3:28-41, 1955. 

Andreas Prohl. Projection and quasi- compressibility methods for solving the incompressible Navier-Stokes equa- 
tions. Advances in Numerical Mathematics. B. G. Teubner, Stuttgart, 1997. 

R. Rannacher. On Chorin's projection method for the incompressible Navier-Stokes equations. In The Navier- 
Stokes Equations II — Theory and Numerical Methods (Oberwolfach, 1991), volume 1530 of Lecture Notes in 
Math., pages 167—183. Springer, Berlin, Germany, 1992. 

J. Shen. On error estimates of some higher order projection and penalty-projection methods for Navier-Stokes 
equations. Numer. Math., 62(l):49-73, 1992. 

J. Shen. On error estimates of projection methods for the Navier-Stokes equations: second-order schemes. Math. 
Comp., 65(215):1039-1065, 1996. 

R. Temam. Sur l'approximation de la solution des equations de Navier-Stokes par la methode des pas fraction- 
naircs ii. Arch. Rat. Mech. Anal, 33:377-385, 1969. 

Roger Temam. Navier-Stokes equations. AMS Chelsea Publishing, Providence, RI, 2001. Theory and numerical 
analysis, Reprint of the 1984 edition. 

L.J. P. Timmermans, P.D. Minev, and F.N. van de Vosse. An approximate projection scheme for incompressible 
flow using spectral elements. Int. J. Numer. Methods Fluids, 22:673-688, 1996. 

Department of Mathematics, Texas A&M University 3368 TAMU, College Station, TX 77843-3368, 
USA. On leave from CNRS, France. 

E-mail address: guermondOmath.tamu.edu 

2 Department of Mathematical and Statistical Sciences, University of Alberta, Edmonton, Alberta 
Canada T6G 2G1. 

E-mail address: minevOualberta. ca 

3 Department of Mathematics, University of Maryland, College Park, MD 20742, USA. 
E-mail address: abnersg<Smath.umd.edu 



[15 

[is; 
[17; 
[is; 
[19; 

[20; 

pi; 
[22; 
[2.3; 
[24; 

[25" 



